%% it computes firm policy of the analytical case used as a guess
u0_l      = (1/nu-chi/A)/(r+delta-mu); %if u0_l<0,error('u0'),end
u1_l      = -(coupon+rho+varphi*delta)/(r+rho+delta);
barphi    = -u1_l;
gamma     =( (mu+rho-0.5*sigma^2)+sqrt((mu+rho-0.5*sigma^2)^2 + 2*sigma^2*(r+delta+rho)) )/sigma^2;

u2_l      = (-(1-phi*a)*u1_l-(1-phi)*varphi)/((1+gamma)*(1-phi*a^(1+gamma)));
b_bar     = u0_l*(1+1/gamma)*(1-phi)/(barphi*(1-phi*a-(1-phi)*varphi));

b_vec      = 0:(b_bar-0)/50:b_bar;
u_second_l = (u2_l*(1+gamma)*gamma*(b_vec/b_bar).^gamma)./b_vec;
u_prime_l  = u1_l+u2_l*(1+gamma)*(b_vec/b_bar).^gamma;
ell_l      = (-u_prime_l./u_second_l)*(dr);
v_l        = u0_l +  u1_l*b_vec + u2_l*(b_vec./b_bar).^gamma.*b_vec;

check=(r+delta-mu)*v_l-(1 -(coupon+rho)*b_vec - (rho+mu)*u_prime_l.*b_vec + sigma^2/2*u_second_l.*b_vec.^2);
if abs(check)>1/100000,error('analytical solution not working'),end

